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ABSTRACT 



We systematically derived hydrodynamic equations and transport coefficients for a class of multi-speed lattice 
Boltzmann models in D dimensions, using the multi-scale technique. The constitutive relation of physical fluid is 
recovered by a modified equilibrium distribution in Maxwell-Boltzmann type. With the use of the rest particles 
and the particle reservoir, we were able to add one degree of freedom into the sound speed of the modeled fluid. 
When the sound speed is tuned small enough, the compressible region of fluid flow can be reached. An example 
2-D model is presented, together with the numerical verification for its transport coefficients. 
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INTRODUCTION 



Is it is well known, the lattice Boltzmann models for 
p^uid dynamics are extremely suitable for the modern 
^^super-computation. The realization is due to some 
Fundamental characteristics of the discrete microscopic 
model: the intrinsic stability, the easiness in dealing 
C^Vith complicated geometric boundaries and the fast 
Cstnd local operations fit for massively parallel process- 
Qng. 



However, the conventional lattice Boltzmann 



Q^ome drawbacks. Apart from the doubt about the effi- 
"■-.ciency of the model on the serial computer frame, the 
~2nain problems result from the unphysical terms ap- 
Lj&earing in the macro-dynamic equations. Specifically, 
^Se modeled fluid lacks of Galilean invariance(presence 
Hsif | a so-called g(p) density dependent factor in the non- 
r+inear advection term) and contains a velocity depen- 
dent part in the pressure term. 
O 

These difficulties restrict the range within which 
the model may be applied to the simulation of fluid 
flow. Accordingly, the incompressible limit is caused 
by two factors. First, the lost Galilean invariance has 
to be recovered by rescaling the time with the coef- 
ficient g(p). This manipulation requires the density 
of the modeled fluid to be set to a constant p . Sec- 
ondly, the sound speed of the modeled fluid is gen- 
erally determined by the lattice geometry in the mi- 
croscopic world. As the velocity vectors of particles 
are identical to the unit base vectors of the regular 
polygon which tesselates the lattice space, the sound 
speed of the model turns out to be order 1. Whereas 
the equilibrium particle distribution function must be 
expanded in the small limit of flow velocity, that is 
u <C c s in the hydrodynamic derivation of the LBE 
model. The Mach number thus has to be kept small, 
which means that the region where the compressibility 



becomes important cannot be reached easily. 

Recently a lattice BGK model was proposed 
(Qian et al. , 1992) to recover the full Navier-Stokes 
equations for the isothermal fluid. The name comes 
from the adoption of the collision operator which had 
been employed in the kinetic BGK equation(Bhatna- 
gar et al., 1954). The operator takes a relaxation form 
in which the particle distribution approaches to the 
equilibrium value in a rate of 1/r. Furthermore, all 
those unphysical terms were eliminated by employing 
the Maxwell-Boltzmann type of equilibrium distribu- 
tion function. Unfortunately the sound speed of the 
model still remained frozen so that it was difficult to 
increase the Mach number. It was claimed that Mach 
number must be less than y/5/2, through a stability 
analysis on the positivity requirement of the particle 
density. More severe conditions would be expected if 
this analysis were done on the spatial gradient of hy- 
drodynamic quantities. Therefore some modification 
of the model is needed to solve the difficulty. 

In this paper, we employ a particle reservoir 
to ensure the mass conservation while the unphysical 
terms are removed by a modified Maxwell-Boltzmann 
type of equilibrium distribution. With such opera- 
tions, we set freely certain member of the t p 's, which 
are factors of the density allocation on different sub- 
lattices. For example, in a 2-D 3-speed model t\ and 
f 2 (allocating factors for speed 1 and 2 particles) are 
designed to recover the isotropy for the fourth order 
tensor, while ^(allocating factor for rest particles) is 
specified as the controller of sound speed of the mod- 
eled fluid. 

We derive the hydrodynamic equations and the 
transport coefficients for a general D-dimensional 
model in Section 2. An example of 2-D model which 
consists of 9 types of particles will be introduced in 
Section 3. Results of numerical calculation of this 



model are presented in Section 4. The conclusion and 
perspective of the model will be discussed in Section 
5. 

2 LATTICE HYDRODYNAMICS 

2.1 Lattice Field And Symmetric Properties 

The lattice field defined in this model is composed of 
several sub-lattices whose base vectors are permuta- 
tions of (±1, . . . , ±1, 0, . . . , 0) in D dimensions. We 
use p, namely, the number of non-zero components in 
the base vector, as the index for different sub-lattices. 
The number of base vectors of each sub-lattice is de- 
noted by bp and can be calculated by 
2PD 

p ~ p\(D-p)\ ■ (1) 
Cpi a (0 <p<D,l<i<b p ,l<a< D) is used to 
denote the a component of possible velocity vectors for 
particles in cell i of sub-lattice p. As c pj 's are identical 
to the base vectors, the following relationship should 
hold, 

(2) 



E 



2 _ 2 _ 
C pia ~ C p — P ■ 



Furthermore, the velocity moments can be explicitly 
given by using the symmetric properties of the regular 
lattice. We have, 
VMP 1 c pia = 0. 



b p c p s: 



VMP 2 J2i C piaCpif3 = —^0 a p. 

VMP 3 Y,i CpiaCpipCpis = 0. 

V]VIP 4 Cpi a CpipCpifiCpiy — ^Pp^a^jS 

Note that the fourth velocity moment will be denoted 
as Tp a p~[S in this paper and ip p and ip p are calculated 
by(Qian, 1990) 

_ ( D + 2 - 3p)b P 4 _ (p-l)b p c$ 

Vp ~ p D (D - 1) ' ^ ~ pD(D - 1) • ( ' 
These expressions can easily be verified by the value 
of Tpim and TJ,i2i2- We define ip p and ip p for rest 
particles as 

t^o = , = . (4) 

2.2 Definition Of Macroscopic Quantities 

The particle distribution in cell i of sub-lattice p at 
time t and lattice site x is designated as N p i(x,t). 
At this stage we introduce a particle reservoir, by 
which the mass conservation may be ensured in the 
microscopic process(Chen et al, 1992). The reser- 
voir works in background at every lattice site and 
delivers products at every time step. The products, 
written as R(x,i) can be either positive or negative. 
When R > 0, they are taken as normal particles which 
may compensate the mass. Otherwise, R is called as 
"particle-killer" which diminish the surplus of parti- 
cles at that site. Hence the macroscopic density and 
momentum flux are defined as, 
DEF 1 p(x,t) = J2 Pii N pi (x,t) + R. 

DEF 2 p(x,t)u(x,t) = Ep,iNpi(x,t)c pi . 



The spatially uniform equilibrium can be reached in 
the lattice field when the macroscopic flow vanishes. 
In this case, the averaged particle density d p at each 
cell of the p sub-lattice satisfies, 

p. (5) 



E 



bp dp 



An allocating factor, which is defined as d p = t p p, is 
used to rewrite equation (5) as 

Y J bpt v = l. (6) 

p 

Obviously b p t p decides the percentage by which the 
particle density on the p sub-lattice contributes to the 
macroscopic mass. We shall see later that t p 's play 
important roles on recovery of isotropy of the fourth 
order tensor and lead into an adjustable sound speed 
for the modeled fluid. In order to ensure the isotropy 
of Tpa/ijs, we must have 



. 



(7) 



We may simply take this requirement as an assumption 
at the present moment. 

2.3 Equilibrium Distribution 

We modify the extended Maxwell-Boltzmann type of 
equilibrium distribution in the following way so that 
both the Galilean invariance and the physical pressure 
term could be achieved in the sequel derivations. We 
write, when u <C 1, 



N, 



(o) 



m +dp[l+j3 (cp ia u a )+hiu 2 + ]^j3' (cpi O! u O! )' 2 ] ,(8) 



for moving particles, 

for rest particles and 
D 

R 



ju 2 - y^bpm 



(9) 



(10) 



p=i 

for the products of the particle reservoir. The 
l3o,hi,/3 and 7 appearing above are given as follows 

A, = ^ D L „ , (11) 



hi 

7 = 



J2 P tpbpCp 
1 

Y^ptpVp ' 

1 - -*oME P tpfp 

-p 



(12) 
(13) 

(14) 



2 A] J2 P i p i fp 

These expressions are derived with the use of VMP 1 
~ VMP 4 and DEF 1 ~ DEF 2. mq is an arbitrary 
constant, introduced to avoid the particle distribution 
to be negative. 

2.4 LBE And Conservation Laws 

The lattice Boltzmann equation is expressed in the en- 
semble averaged form of the micro-dynamic equation 



for particles moving on sub-lattices. It is given as fol- 
lows, 

N pi (x + cpi,t + 1) - N pi (x,t) = A pi (N) , (15) 
with the Boltzmann approximation applied to the col- 
lision term Ap 8 -. Thus the LBE exactly describes the 
two-step motion of particles in the micro-world. The 
first step is the streaming period during which the local 
particle distribution is transported to one of its near- 
est neighbors. The second one is the collision period 
during which the particle distribution relaxes to its 
equilibrium value which locally depends on the macro- 
scopic density and momentum. A pj - may be linearized 
in the form 

A pi (N) = A pi (N^) + A%(N 9j - < 0) ) . (16) 
From the definition of the collision process, we con- 
clude that the first term will vanish as is already 
the equilibrium value. A p q 'j is the so-called collision 
operator which gives the probability for particle tran- 
sition from the cell j of sub-lattice q to the cell i of 
sub-lattice p. The operator of BGK type takes the 
simplest but most efficient form as 

A^=-U pq 8 lJ . (17) 

The mass and momentum are conserved globally in 
the streaming period while they are kept constant lo- 
cally in the collision process. Such conservation rela- 
tion may be described as follows, 



2.5.1 On e Order 



A pi (N) = 0,J2 A pi( N ) c pi 



. 



(18) 



p,i p,i 

These relations would lead to the macroscopic conser- 
vation laws by substituting the equilibrium distribu- 
tion function and rescaling time and space properly. 

2.5 Multi-Scale Technique 

In the long wavelength and low frequency limit, the 
LBE can be expanded as follows, 
1 „ 

dtN p i + c pia d a N pi + -d t N pi + c pia d t d a N pi 
1 Z 

+ -c pia c p if j d a df j N pi = A p i(N) . (19) 

The multi-scale technique is applied in the following 
way. First we find a proper quantity e<l and write 
N p i perturbatively with the use of e, 



N, 



#i 0) + cNP 



Note that the perturbative part must satisfy the fol- 
lowing constraints, 

E< 1) = ' E< 1)c - = ' ( 21 ) 

p,i p,i 

due to the definition of the mass and momentum and 
the conservation relation (18). The time and space 
differential operator will also be replaced by 

d t ^(ed tl + e 2 d t2 ),d a ^ed a . (22) 
When the expressions stated above are substituted 
into equation (19), we shall pick up terms on the same 
order and formalize them in different equations. 



As terms on this ord 
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dnN^ + c pia d a N 
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(23) 



Next the equation (23) and its velocity moment will be 
integrated in the discrete velocity space. Employing 
VMP 1 ~ VMP 4, DEF 1 ~ DEF 2 and equation 
(18), (21), two macroscopic equations appear below, 

dtiP+d p (pup) = 0, (24) 

dn(pu a ) + dp{pu a up) = -^-d a p . (25) 

Po 

2.5.2 On e 2 Order 

The terms kept on this order are formalized as 
follows , 

d n N^ + c pia d a N^ + d t2 N^ + \^N^ + 

CpiadtldaN^ + -CpiaCpipdadpN^ = . (26) 

The integration would be carried out in the same way 
as the first order equation. The integrated terms must 
be identified with the relations mentioned in the pre- 
vious section, together with equation (24) and (25) 

themselves. Furthermore, the expression for N ( p { } can 
be found from equation (23) and the order would be 
kept only to e 2 u in the reduced equations. They are 
described below, 

d t 2P = , (27) 
d t 2(pu a ) + (t - -)-3-<9 a [d 7 (pu 7 )] - (t- -) y^tp^pp 

/3o5 /3 [(9 7 (/9M 7 )(5„ / 3 + d a (puf}) + dp(pu a )] = . (28) 
2.6 Macroscopic Equations 

When equations (24), (25) and (27), (28) are re- 
combined by using the relation (22), the full Navier- 
Stokes equations emerge as follows, 
d t p + dp(pu p ) = 0, (29) 

dt{pu a ) + dp(pu a up) = 

-c]d a p + r]dj(pu a ) + C<9a[3 7 (/ou 7 )] . (30) 

the transport coefficients are given by 

1 (31) 



p 

c = (r - bi^o^p^p - h ■ 



(32) 



(33) 



The derivation of hydrodynamic equations is 
completed at this point. However, the allocating 
factors(tp's) are still to be determined by relation (6) 
and (7) so that the transport coefficients can be cal- 
culated explicitly. In the next section we shall give a 
2 dimensional model and show r and t p 's influence on 
the viscosity and sound speed of the modeled fluid. 







1 y 






C22 


C12 




C21 






/coi 


X 


C13 






Cll 




C25 


C14 


C24 



Figure 1: Lattice geometry for a 3-speed model in two 
dimensions. The arrows represent the velocity vectors 
of different types of particles. 



3 AN EXAMPLE 2-D MODEL 

We present a 2-D model which at most contains 9 kinds 
of particles at one lattice site. The lattice employed 
in the model is sketched in figure 1. We may imme- 
diately obtain the set of parameters concerning with 
the lattice geometry. They are classified by the lattice 
index p, whose possible values would be 0,1,2 in this 
case, 

r & = i r ^o = o r ^o = o 

I 6i=4 , * = 2 , < <Pi = . (34) 
[ & 2 = 4 [ i> 2 = -8 { ip 2 = 4 
In order to calculate the equilibrium distribution and 
the transport coefficients, t p 's have to be decided, as 
mentioned before, with equation (6) and (7). The so- 
lution is 

'Ii-^o • ( 35 ) 
~ ~ 2 °~ 

From these expressions, we first give the formulas for 
the calculation of equilibrium distribution of different 
particles and the products of the particle reservoir as 
follows , 

(o) 
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N. 
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t p, (36) 

2 ] ,(37) 

Vl ,(38) 



N\"> = m + di+ p[-(u a c lia )+ -(u a c lia )- - -u 

6 1 6 

(0) 



m +d2+p[-^-(u a C2ia)+l(UaC2ia) 2 ~ 
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R. 



-8m - -pu 2 . (39) 

The average density of the speed-1 particle is kept to 
be 4 times of that of the speed-2 particle so as to result 
in an isotropic fourth order tensor. Referring to those 
macroscopic transport coefficients, they can also be 
written explicitly as functions of to and r, 



C = (r-|)(|-^). 



(41) 
(42) 



r(l - to) 



(40) 



The sound speed can be tuned to a suitable value by 
properly selecting the t . We observed that the de- 
pendence on (1 — to) borne by c s in equation (40) 
is physically meaningful. When t = 0, or there are 
no rest particles on the lattice field, the sound speed 
would reach its maximum value. On the contrary, as 
t = 1, the expression indicates a zero sound speed. 
Because in the latter case, variations in the density 
would be absorbed completely into the rest particles 
so that no sound wave could be transported. The bulk 
viscosity's dependence on c s is nontrivial, as indicated 
in the analysis of the standard fluid mechanics(Landu 
and Lifshitz, 1959). Finally the positivity of the shear 
viscosity requires that r should always be larger than 
l 

2 • 

4 NUMERICAL EXPERIMENTS 

4.1 Measurement Of Transport Coefficients 

The method used to measure the transport coefficient 
of the modeled fluid is described as follows. First a 
periodically perturbative field of mass and momentum 
should be set up artificially, then the relaxing process 
would be measured and spatially Fourier-transformed. 
With the use of the curve fit technique, the relation 
between the wave vector k and the angular frequency 
ijj could be found out, from which the values of various 
transport coefficients can be generalized. 

For simplicity, the wave vector k is set to be 
(k, 0) while the uniform motion of the modeled fluid is 
characterized by (uo, 0). Thus the perturbated density 
and momentum can be written as 
p=p + 6pe i (- ut+k ^ , (43) 

pu x = pu + 8j x e^ t+k ^ , (44) 

p Uy = 6j y e^ t+k ^ . (45) 

Here bp, 8j x and 6j y are small quantities which de- 
cide the amplitude of the perturbation. When these 
expressions are substituted into the macro-dynamic 
equations and higher order terms are dropped, a set of 
linearized hydrodynamic equations emerge as follows, 
iu> bp + ik 6j x = , (46) 

iuj 8j x + ik (u 5p + 28j x )gu = -ik c 2 s 8p- 

k 2 ^ + C)6j x , (47) 

iui 8j y + ikgu 8j y = -k 2 r]8j y . (48) 
We have deliberately left an additional coefficient be- 
fore the advection term. If the Galilean invariance is 
recovered by the model, it should be 1. In the sequel, 
we shall use proper combination of equations (46), (47) 
and (48) to measure the sound speed c s , the shear vis- 
cosity rj, the bulk viscosity £ and the advection coeffi- 
cient g of the modeled fluid. 

4.1.1 Sound Speed 

To measure the sound speed, we need to set u = 
and neglect the effect of viscosity. We use equations 



(46) and (47) to obtain the following expression if only 
the density perturbation is considered, 



Here u> p denote the angular frequency for the damping 
of the density perturbations. The numerical result is 
shown in figure 2. We see that the difference to the 
theoretical prediction is trivial. 
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Figure 2: Sound speed of the modeled fluid as a func- 
tion of to- The solid line is the theoretical prediction. 
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Figure 3: Shear viscosity of the modeled fluid as a 
function of r. The solid line is the theoretical predic- 
tion. 



approaching to the equilibrium state through the col- 
lision process. When the Knudsen number increases, 
the mean free path of particles would get to be compa- 
rable with the system size so that the hydrodynamic 
mode would be broken down eventually. 



4.1.2 Shear and Bulk Viscosity 

The expression of the shear viscosity can easily 
be obtained from equation (48) with u = 0, 



k 2 



(50) 



Here uij y refers to the angular frequency for the damp- 
ing of the momentum perturbation in the ^-direction. 
In the same spirit we can derive the formula for bulk 
viscosity from equation (46) and (47), by considering 
the damping of the momentum perturbation in the IE- 
direction, 



k 2 



T]. 



(51) 



Here we assume that the sound speed and the shear 
viscosity have already been known. The comparison 
between results of the numerical measurements and 
the theoretical predictions is shown in figure 3 and 
figure 4, respectively. 

There obviously exists deviation from theoreti- 
cal solution in the shear viscosity when the relaxation 
time becomes large enough. This is due to the large 
Knudsen number which is characterized by a slower 
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Figure 4: Bulk viscosity of the modeled fluid as a func- 
tion of sound speed. The solid line is the theoretical 
prediction. 
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5 CONCLUSION 

We extended the range of application of a multi-speed 
LBE model by introducing one degree of freedom in 
the sound speed of the modeled fluid. The analytical 
results seem to be physically meaningful while the nu- 
merical verification show good agreement with them. 

When the sound speed of the modeled fluid is 
tuned to be comparable with the flow velocity, the 
Mach number could be much higer than the conven- 
tional scheme. This would make it easier for the model 
to simulate the supersonic flows. Moreover, the con- 
troller of the sound speed, namely to, need not to be 
constant. We may design it to be properly depen- 
dent on thermodynamic variables(Chen et al, 1993) 
so that some realistic isentropic flow problems can be 
simulated either. 

On the other hand, because the temperature has 
not been defined in the model, many practical appli- 
cations are not amenable so far. We think that the 
difficulty would be overcome in the near future. 



4.1.3 Advection Coefficient 

We initialize the flow motion with a nonvanishing 
u , and switch off the viscous effects in the mean time. 
Perturbations are added to the momentum density in 
the ^-direction. From equation (48) we obtain, 

^-ib^ , (52) 

The result is printed in figure 5, which can be taken 
as a complete proof of the recovery of the Galilean 
invariance in the modeled fluid. 




